##########################################################################################

library(data.table)
library(optparse)
library(ggplot2)
library("FactoMineR")
library("factoextra")

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--sample_list_public_file"), type = "character"),
    make_option(c("--raw_tcga_count_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--gtf_for_len_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    combine_tmm_file <- "~/20220915_gastric_multiple/dna_combinePublic/mRNA/CombineCounts.TCGA_NJMU.FilterLowExpression.TMM.tsv"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/mRNA"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
rsem_file <- opt$rsem_file
sample_list_public_file <- opt$sample_list_public_file
raw_tcga_count_file <- opt$raw_tcga_count_file
out_path <- opt$out_path
gtf_file <- opt$gtf_file
gtf_for_len_file <- opt$gtf_for_len_file

dir.create(out_path , recursive = T)

##########################################################################################

normalized_counts <- data.frame(fread(combine_tmm_file))
colnames(normalized_counts) <- gsub( "[.]" , "-" , colnames(normalized_counts) )

##########################################################################################

dat <- normalized_counts
njmu_sample <- grep( "TCGA" , colnames(dat) , invert = T , value  = T)
tcga_sample <- grep( "TCGA" , colnames(dat) , value  = T)

disease <- sapply(strsplit(sapply(strsplit(colnames(dat)[-ncol(dat)] , "_") , "[" , 2) , "-"  ) , "[" , 1)
batch <- batch <- ifelse( colnames(dat) %in% njmu_sample , "NJMU"  , "TCGA")[-ncol(dat)]

##########################################################################################
## 提取Normal的离群样本
id <- "Normal"
tmp <- dat[,which(disease==id)]
from <- batch[which(disease==id)]
df.pca <- PCA(t(tmp), graph = FALSE)

## 基于PCA图观察离群点
## CombineCounts.TCGA_NJMU.FilterLowExpression.TMM.PCA_Batch_Normal.pdf
normal_rm_sample <- names(which(df.pca$ind$coord[,2] > 60))

id <- "IGC"
tmp <- dat[,which(disease==id)]
from <- batch[which(disease==id)]
df.pca <- PCA(t(tmp), graph = FALSE)
igc_rm_sample <- names(which(df.pca$ind$coord[,2] > 100))

id <- "DGC"
tmp <- dat[,which(disease==id)]
from <- batch[which(disease==id)]
df.pca <- PCA(t(tmp), graph = FALSE)
dgc_rm_sample <- names(which(df.pca$ind$coord[,2] > 100))

##########################################################################################
## 去除离群样本
dat_qc <- dat[,!colnames(dat) %in% c(normal_rm_sample , igc_rm_sample , dgc_rm_sample)]

dat <- dat_qc
njmu_sample <- grep( "TCGA" , colnames(dat) , invert = T , value  = T)
tcga_sample <- grep( "TCGA" , colnames(dat) , value  = T)

disease <- sapply(strsplit(sapply(strsplit(colnames(dat)[-ncol(dat)] , "_") , "[" , 2) , "-"  ) , "[" , 1)
batch <- batch <- ifelse( colnames(dat) %in% njmu_sample , "NJMU"  , "TCGA")[-ncol(dat)]

pca.plot = function(dat , col , out_file){

    df.pca <- PCA(t(dat), graph = FALSE)
    fviz_pca_ind(df.pca,
       geom.ind = "point",
       col.ind = col ,
       addEllipses = TRUE,
       legend.title = "Groups"
    )
}

## 现在的批次
out_file <- paste0( out_path , "/CombineCounts.TCGA_NJMU.FilterLowExpression.TMM.PCA_Batch.QC.pdf" )
pdf(out_file)
pca.plot( dat[,-ncol(dat)], factor(batch) , out_file)
dev.off()

id <- "Normal"
out_file <- paste0( out_path , "/CombineCounts.TCGA_NJMU.FilterLowExpression.TMM.PCA_Batch_",id,".QC.pdf" )
pdf(out_file)
tmp <- dat[,which(disease==id)]
from <- batch[which(disease==id)]
## 现在的批次
pca.plot( tmp, factor(from) , out_file)
dev.off()

id <- "IGC"
out_file <- paste0( out_path , "/CombineCounts.TCGA_NJMU.FilterLowExpression.TMM.PCA_Batch_",id,".QC.pdf" )
pdf(out_file)
tmp <- dat[,which(disease==id)]
from <- batch[which(disease==id)]
## 现在的批次
pca.plot( tmp, factor(from) , out_file)
dev.off()

id <- "DGC"
out_file <- paste0( out_path , "/CombineCounts.TCGA_NJMU.FilterLowExpression.TMM.PCA_Batch_",id,".QC.pdf" )
pdf(out_file)
tmp <- dat[,which(disease==id)]
from <- batch[which(disease==id)]
## 现在的批次
pca.plot( tmp, factor(from) , out_file)
dev.off()
